The invention relates to the calibration of optical techniques for imaging and spectroscopy of, e.g., biological tissue.
Recently there has been significant interest in using optical radiation for imaging within highly scattering media, such as biological tissue. Photons travel within the highly scattering media along a distribution of paths, of which very few are straight. Thus, directing light into the highly scattering media and subsequently detecting the diffuse light emitted from the media provides information about local variations in the scattering and absorption coefficients. Such information can identify, for example, an early breast or brain tumor, a small amount of bleeding in the brain, or an early aneurysm. Furthermore, for example, multiple wavelengths can be used to spectroscopically determine local tissue concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin (Hb) in tissue, which may be in response to some stimulus, e.g., a drug. For a general description of such applications, see, e.g., A. Yodh and B. Chance in xe2x80x9cSpectroscopy and Imaging with Diffusing Light,xe2x80x9d Physics Today, pp. 34-40 (March 1995).
If the spatially varying optical properties of the highly scattering media are known, photon propagation within the media can be calculated numerically. The numerical calculation is simplified when scattering is much larger than absorption, in which case the photon propagation can be approximated as a diffusion process. This condition is typically satisfied in biological tissue in the spectral range of about 700 nm to 900 nm. The numerical calculation gives the distribution of light inside the tissue, and is usually referred to as the xe2x80x9cforward calculation.xe2x80x9d For a sample being imaged, however, the xe2x80x9cinverse calculationxe2x80x9d needs to be solved, i.e., deducing the sample""s optical properties from the diffuse light measurements. Numerical techniques for performing the inversion include singular value decomposition (SVD), simultaneous iterative reconstruction technique (SIRT), K-space diffraction tomography, and using an extended Kaman filter. For a general review of techniques for the forward and inverse calculations, see, e.g., S. R. Arridge in xe2x80x9cOptical tomography in medical imaging,xe2x80x9d Inverse Problems, 15:R41-R93 (1999).
In diffuse optical tomography (DOT), multiple sources sequentially direct light into tissue at spatially separated locations, and for each such source, multiple detectors on the tissue measure the diffuse light emitted from the sample at spatially separated locations. For every source-detector pair, one measures the local transmittance of the diffuse light, i.e., the ratio of the diffuse radiance measured by the detector and the incident radiance from the source. The measurements provide the input information for the inverse calculation. However, the measurements can include various uncertainties caused by, for example, fluctuations in the source power, variations in the detector gain, and coupling variations at the source-tissue interface as well as the tissue-detector interface.
To minimize the uncertainties, DOT systems are typically calibrated with initial measurements for a known sample, and the calibration is used to correct subsequent measurements for imaging an unknown sample. Unfortunately, coupling at the source-tissue interface and the tissue-detector interface can vary from measurement to measurement because of, for example, movement or perspiration of the patient, or movement of an optical fiber that forms part of a source or detector. Thus, the results of the inverse calculation can include systematic errors caused by measurement variations that are independent of the tissue optical properties of interest.
The systematic errors can also limit absolute spectroscopic measurements of optical properties at a particular spatial location, e.g., the absolute, rather than relative, values of absorption and scattering.
The invention features a calibration method for diffuse optical measurements that corrects transmittance measurements between a source land a detector for factors unrelated to sample properties. For imaging applications, the corrected transmittance measurements can be subject to an inverse calculation to determine spatial variations in the optical properties of the sample, i.e., to xe2x80x9cimagexe2x80x9d the sample. For spectroscopic applications, the corrected transmittance measurements can be used to determine absolute values for the optical properties of the sample in a particular spatial region at multiple wavelengths, e.g., to determine the absolute concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin (Hb). The calibration method is based on the same set of transmittance measurements that are subsequently corrected by the calibration and used in imaging and/or spectroscopy applications. The accuracy of the subsequent results is thus not subject to uncertainties caused by a delay between calibration and sample measurements.
The calibration method involves a forward calculation for each of multiple source-detector pairs based on an approximate model of the sample, and a minimization of an expression that depends on the forward calculations and the transmittance measurements to determine self-consistent coupling coefficients for every source-detector pair. Once the coupling coefficients have been determined, they can be used to correct the transmittance measurements. If desired, an inverse calculation can be performed on the corrected sample measurements to determine spatial variations in the optical properties of the sample. If necessary, the calibration can be repeated and iteratively improved, whereby the optical properties determined by the inverse calculation in an earlier iteration are used to improve the sample model for the forward calculation in a subsequent iteration.
In general, in one aspect, the invention features a system for making optical measurements. The system includes at least two optical sources which during operation couple optical radiation into a sample at spatially separated locations and at least two optical detectors positioned to receive optical radiation emitted from the sample at spatially separated locations in response to the optical radiation from the sources, and an analyzer.
The signal g(i,j) produced by the jth detector in response to the optical radiation from the ith source can be expressed as g(i,j)=SiDjf(i,j), where f(i,j) depends only on the properties of the sample, Si is a coupling coefficient for the ith source, and Dj is a coupling coefficient the jth detector. During operation, the analyzer calculates the value of the product SlDk for at least one of the source-detector pairs based on the signals produced by the detectors and simulated values of f(i,j) corresponding to a model of the optical properties of the sample. The sources, for example, can provide continuous-wave radiation, in which case g(i,j), f(i,j), Si, and Dj are all real-valued. Alternatively, the sources can provide modulated CW radiation, or short temporal pulses of optical radiation (e.g., less than about 1 to 100 ps). In these latter cases, g(i,j), f(i,j), Si, and Dj can be complex, or more generally they can be vectors representing multiple values (e.g., transmittance and temporal delay).
Embodiments of the system can also include any of the following features.
The analyzer can calculate the value of the product SiDj for every source-detector pair based on the detector signals and the simulated values of f(i,j). The analyzer can further calculate experimental values of f(i,j) based on the calculated values of SiDj and the signals g(i,j) using the expression g(i,j)=SiDjf(i,j), and then perform an inverse calculation on the experimental values for f(i,j) to determine spatial variations in at least one optical property of the sample. The optical property or properties can be the absorption coefficient, the reduced scattering coefficient, or both. The analyzer can modify the model of the sample based on the determined spatial variations, and then repeat the calculation of the values of the product SiDj for every source-detector pair using the modified model. The initial model can correspond to the sample being homogeneous.
The analyzer can simulate values of f(i,j) according to the expression:       f    ⁡          (              i        ,        j            )        ∝                    3        ⁢                  μ          s          xe2x80x2                            4        ⁢        π        ⁢                  "LeftBracketingBar"                      r            ij                    "RightBracketingBar"                      ⁢          xe2x80x83        ⁢          exp      ⁡              [                              -                                          (                                  3                  ⁢                                      μ                    s                    xe2x80x2                                    ⁢                                      μ                    a                                                  )                                            1                2                                              ⁢                      "LeftBracketingBar"                          r              ij                        "RightBracketingBar"                          ]            
where |rij| is the distance between the ith source and the jth detector. This expression corresponds to an infinite, homogenous medium.
The analyzer can calculate the value of the product SlDk by minimizing the expression:             F      ⁡              (                              S            l                    ⁢                      D            k                          )              =                  ∑                  i          =          1                          N          s                    ⁢                        ∑                      j            =            1                                N            d                          ⁢                              (                                                                                L                    ⁡                                          (                                              i                        ,                        j                        ,                        k                        ,                        l                                            )                                                                                                  S                      l                                        ⁢                                          D                      k                                                                      ⁢                                  xe2x80x83                                ⁢                                  f                  ⁡                                      (                                          i                      ,                      j                                        )                                                              -                              g                ⁡                                  (                                      i                    ,                    j                                    )                                                      )                    2                      ,
where
L(i,j,k,l)=AsikAdjl
            A      s      ik        =                            N          s                          N          d                    ⁢                        ∑                      j            =            1                                N            d                          ⁢                              g            ⁡                          (                              i                ,                j                            )                                                          f              ⁡                              (                                  i                  ,                  j                                )                                      ·                                          ∑                                  ii                  =                  1                                                  N                  s                                            ⁢                                                                    g                    ⁡                                          (                                              ii                        ,                        j                                            )                                                                            f                    ⁡                                          (                                              ii                        ,                        j                                            )                                                                      ·                                                      f                    ⁡                                          (                                              ii                        ,                        k                                            )                                                                            g                    ⁡                                          (                                              ii                        ,                        k                                            )                                                                                                                          A      d      jl        =                            N          d                          N          s                    ⁢                        ∑                      i            =            1                                N            s                          ⁢                              g            ⁡                          (                              i                ,                j                            )                                                          f              ⁡                              (                                  i                  ,                  j                                )                                      ·                                          ∑                                  jj                  =                  1                                                  N                  d                                            ⁢                                                                    g                    ⁡                                          (                                              i                        ,                        jj                                            )                                                                            f                    ⁡                                          (                                              i                        ,                        jj                                            )                                                                      ·                                                      f                    ⁡                                          (                                              l                        ,                        jj                                            )                                                                            g                    ⁡                                          (                                              l                        ,                        jj                                            )                                                                                                              
and NS is the number of sources and ND is the number of detectors.
Furthermore, when the model corresponds to the sample being homogeneous, the analyzer can calculate at least one of the absorption coefficient xcexca and the reduced scattering coefficient xcexcxe2x80x2s by minimizing F(SlDk) with respect lo the product SlDk and the at least one of the absorption and scattering coefficients. F(SlDk) implicitly depends on xcexca and xcexcxe2x80x2s through f(i,j). Similarly, the analyzer can calculate both of the absorption and scattering coefficients by minimizing F(SlDk) with respect to the product SlDk and the absorption and scattering coefficients.
The analyzer can also calculate the product SmDn for every source-detector pair by minimizing F(SmDn). Alternatively, the analyzer can calculates the product SiDj for every remaining source-detector pair according to             S      i        ⁢          D      j        =                    L        ⁡                  (                      i            ,            j            ,            k            ,            l                    )                                      S          l                ⁢                  D          k                      .  
In another aspect, the invention features a method for calibrating an optical measurement system. The optical measurement system includes at least two optical sources and at least two optical detectors. The sources couple optical radiation into a sample at spatially separated locations and the detectors are positioned to receive optical radiation emitted from the sample at spatially separated locations and generate signals in response to the optical radiation from the sources.
The methods includes the steps of: providing the signals generated by the detectors, wherein the signal g(i,j) generated by the jth detector in response to the optical radiation from the ith source can be expressed as g(i,j)=SiDjf(i,j), where f(i,j) depends only on the properties of the sample, Si is a coupling coefficient for the ith source, and Dj is a coupling coefficient for the jth detector; and calculating the value of the product SlDk for at least one of the source-detector pairs based on the signals generated by the detectors and simulated values of f(i,j) corresponding to a model of the optical properties of the sample.
The sources in the optical measurement system, for example, can provide continuous-wave radiation, in which case g(i,j), f(i,j), Si, and Dj are all real-valued. Alternatively, the sources can provide modulated CW radiation, or short temporal pulses of optical radiation (e.g., less than about 1 to 100 ps). In these latter cases, g(i,j), f(i,j), Si, and Dj can be complex, or more generally they can be vectors representing multiple values (e.g., transmittance and temporal delay).
Embodiments of the calibration method can also include any of the following features.
The calibration method can calculate the value of the product SiDj for every source-detector pair based on the detector signals and the simulated values of f(i,j). The method can further include calculating experimental values of f(i,j) based on the calculated values of SiDj and the signals g(i,j) using the expression g(i,j)=SiDjf(i,j), and then performing an inverse calculation on the experimental values for f(i,j) to determine spatial variations in at least one optical property of the sample. The optical property or properties can be the absorption coefficient, the reduced scattering coefficient, or both. Furthermore, the method can include modifying the model of the sample based on the determined spatial variations, and then repeating the calculation of the values of the product SiDj for every source-detector pair using the modified model. The initial model can correspond to the sample being homogeneous.
The simulated values of f(i,j) can be calculated according to the expression:       f    ⁡          (              i        ,        j            )        ∝                    3        ⁢                  μ          s          xe2x80x2                            4        ⁢        π        ⁢                  "LeftBracketingBar"                      r            ij                    "RightBracketingBar"                      ⁢          xe2x80x83        ⁢          exp      ⁡              [                              -                                          (                                  3                  ⁢                                      μ                    s                    xe2x80x2                                    ⁢                                      μ                    a                                                  )                                            1                2                                              ⁢                      "LeftBracketingBar"                          r              ij                        "RightBracketingBar"                          ]            
where |rij| is the distance between the ith source and the jth detector. This expression corresponds to an infinite, homogenous medium.
The method can calculate the value of the product SlDk by minimizing the expression:             F      ⁡              (                              S            l                    ⁢                      D            k                          )              =                  ∑                  i          =          1                          N          s                    ⁢                        ∑                      j            =            1                                N            d                          ⁢                              (                                                                                L                    ⁡                                          (                                              i                        ,                        j                        ,                        k                        ,                        l                                            )                                                                                                  S                      l                                        ⁢                                          D                      k                                                                      ⁢                                  xe2x80x83                                ⁢                                  f                  ⁡                                      (                                          i                      ,                      j                                        )                                                              -                              g                ⁡                                  (                                      i                    ,                    j                                    )                                                      )                    2                      ,
where
L(i,j,k,l)=AsikAdjl
            A      s      ik        =                            N          s                          N          d                    ⁢                        ∑                      j            =            1                                N            d                          ⁢                              g            ⁡                          (                              i                ,                j                            )                                                          f              ⁡                              (                                  i                  ,                  j                                )                                      ·                                          ∑                                  ii                  =                  1                                                  N                  s                                            ⁢                                                                    g                    ⁡                                          (                                              ii                        ,                        j                                            )                                                                            f                    ⁡                                          (                                              ii                        ,                        j                                            )                                                                      ·                                                      f                    ⁡                                          (                                              ii                        ,                        k                                            )                                                                            g                    ⁡                                          (                                              ii                        ,                        k                                            )                                                                                                                          A      d      jl        =                            N          d                          N          s                    ⁢                        ∑                      i            =            1                                N            s                          ⁢                              g            ⁡                          (                              i                ,                j                            )                                                          f              ⁡                              (                                  i                  ,                  j                                )                                      ·                                          ∑                                  jj                  =                  1                                                  N                  d                                            ⁢                                                                    g                    ⁡                                          (                                              i                        ,                        jj                                            )                                                                            f                    ⁡                                          (                                              i                        ,                        jj                                            )                                                                      ·                                                      f                    ⁡                                          (                                              l                        ,                        jj                                            )                                                                            g                    ⁡                                          (                                              l                        ,                        jj                                            )                                                                                                              
and NS is the number of sources and ND is the number of detectors.
Furthermore, when the model corresponds to the sample being homogeneous, the method can include calculating at least one of the absorption coefficient xcexca and the reduced scattering coefficient xcexcxe2x80x2s by minimizing F(SlDk) with respect to the product SlDk and the at least one of the absorption and scattering coefficients. F(SlDk) implicitly depends on xcexca and xcexcxe2x80x2s through f(i,j). Similarly, the method can include calculating both of the absorption and scattering coefficients by minimizing F(SlDk) with respect to the product SlDk and the absorption and scattering coefficients.
The method can also include calculating the product SmDn for every source-detector pair by minimizing F(SmDn). Alternatively, the method can calculates the product SiDj for every remaining source-detector pair according to             S      i        ⁢          D      j        =                    L        ⁡                  (                      i            ,            j            ,            k            ,            l                    )                                      S          l                ⁢                  D          k                      .  
In another aspect, the invention features a computer-readable medium which causes a processor to perform the steps of any of the embodiments of the calibration method described above.
Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. Although methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present invention, the suitable methods and materials are described below. All publications, patent applications, patents, and other references mentioned herein are incorporated by reference in their entirety. In case of conflict, the present specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and not intended to be limiting.
Embodiments of the invention include many advantages. For example, the new calibration methods eliminate errors caused by using separate data sets to calibrate and image samples. Where the sample is human tissue, such errors can result from patient or fiber movement, or patient perspiration. Furthermore, the calibration method obviates the need to minimize fluctuations in the source outputs and detector gains. Diffuse optical measurement systems employing the new calibration methods can therefore provide more accurate and robust results.
Other features and advantages of the invention will be apparent from the following detailed description, and from the claims.